#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<omp.h>
#include<time.h>
#include<ctype.h>
#include<mpi.h>

void comprovacao(int J_ORDER, double x[J_ORDER], int J_ROW_TEST, double **matriz, double b[J_ORDER], int it) {
	int i;
	double valor = 0;
        
    for(i = 0; i < J_ORDER; i++) {
    	valor += matriz[J_ROW_TEST][i]*x[i];
    }
        
    printf("---------------------------------------\n");
    printf("Iterations: %d\n", it);
    printf("RowTest: %d => [%lf] =? %lf\n", J_ROW_TEST, valor, b[J_ROW_TEST]);
    printf("---------------------------------------\n");
        
}

//verificacao de convergencia
int verifica(int J_ORDER, double **aEstrela) {
    int i, j;
    double x[J_ORDER], max;

    for(i = 0; i < J_ORDER; i++) {
        x[i] = 0;
    }

    for(i = 0; i < J_ORDER; i++) {
        for(j = 0; j < J_ORDER; j++) {
            if(i != j){
                x[i] += aEstrela[i][j];
            }
        }
    }
    max = x[0];

    for(i = 1; i < J_ORDER; i++) {
       if(max < x[i])
           max = x[i];
    }

    if(max < 1)
        return 1;

    else
        return 0;
}

int main(int argc, char *argv[]) {
	int J_ORDER, J_ROW_TEST, J_ITE_MAX, i, j;
	double J_ERROR, **matriz, **aEstrela; 
	clock_t inicio = 0, fim = 0;
    inicio = clock();  
  	FILE *pFile;
  	
  	pFile = fopen ( "matriz250.txt" , "rb" );
  	if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  	rewind (pFile); //volta para o início do arquivo
  	
  	fscanf(pFile, "%d", &J_ORDER);
  	fscanf(pFile, "%d", &J_ROW_TEST);
  	fscanf(pFile, "%lf", &J_ERROR);
  	fscanf(pFile, "%d", &J_ITE_MAX);
  	
  	double b[J_ORDER], x0[J_ORDER];
  
  	matriz = malloc( J_ORDER * sizeof(double*) );
	#pragma omp parallel for
  	for(i = 0; i < J_ORDER; i++) {
    	matriz[i] = malloc( J_ORDER * sizeof(double) );
  	}
  	
  	aEstrela = malloc( J_ORDER * sizeof(double*) );
  	#pragma omp parallel for
    for(i = 0; i < J_ORDER; i++) {
        aEstrela[i] = malloc( J_ORDER * sizeof(double) );
    }	
  	
  	//leitura da matriz
  	for(i = 0; i < J_ORDER; i++) {
    	for(j = 0; j < J_ORDER; j++) {
        	fscanf(pFile, "%lf", &matriz[i][j]);
        }
  	}
  	
  	//leitura de b
  	for(i = 0; i < J_ORDER; i++) {
        fscanf(pFile, "%lf", &b[i]);
        x0[i] = 0;
  	}
  	
  	//calculo do A*
  	for(i = 0; i < J_ORDER; i++) {
    	for(j = 0; j < J_ORDER; j++) {
        	if(i != j) {
            	aEstrela[i][j] = (matriz[i][j] / matriz[i][i])*(-1);
        	}
        }
    }

    if(verifica(J_ORDER, aEstrela) != 1) {
        return;
    }
    else {
        double bEstrela[J_ORDER];
        
        //calculo do b*
        for(i = 0; i < J_ORDER; i++) {
            bEstrela[i] = b[i]/matriz[i][i];
        }
        
        int ret, numTasks, myRank;
		MPI_Status stat;
	
		ret = MPI_Init(&argc, &argv);
		if(ret != MPI_SUCCESS)
			exit(1);
		
		MPI_Comm_size(MPI_COMM_WORLD, &numTasks);
		MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
		
		//processo pai
		if(myRank == 0) {
			int continuar = 1, it = 0, tam, tamUltimo;
			double e = 0, max1, max2, x[J_ITE_MAX][J_ORDER];// e = erro relativo\/max1 auxiliar para ver qual a maior diferença entre Xk+1 e Xk e max2 é o maior X
			
			//numero de linhas da matriz A* que sera tratado nos filhos
			tam = J_ORDER/(numTasks-1);
			
			//se a ordem nao for divisivel pelo numero de processos filhos o ultimo filho a ler a matriz ira ler um pedaco menor
			if(J_ORDER%(numTasks-1) != 0) {
				tam = tam + 1;
				tamUltimo = abs(J_ORDER - ((numTasks - 2)*(tam)));	
			}
			else {
				tamUltimo = tam;
			}
			
			for(i = 0; i < J_ORDER; i++){
			    x[0][i] = x0[i];
		    }
		   
		   	//enviando as informacoes para os filhos
		    for(i=1; i<numTasks; i++) {
		    	MPI_Send(&x, J_ITE_MAX*J_ORDER, MPI_DOUBLE, i, 1, MPI_COMM_WORLD);
		    	MPI_Send(&continuar, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
		    	MPI_Send(&it, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
		    	if(i == numTasks-1) {
		    		MPI_Send(&tamUltimo, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
		    	}
		    	else {
		    		MPI_Send(&tam, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
		    	}
		    }
		    
		    while(continuar == 1) {  //enquanto erro é menor...
				for(i = 1; i < numTasks;i++){
					MPI_Recv(&it, 1, MPI_INT, i, 1, MPI_COMM_WORLD, &stat);
					if(i == numTasks-1) {
						MPI_Recv(&x[it][(i-1)*tamUltimo],tamUltimo,MPI_DOUBLE,i,1,MPI_COMM_WORLD,&stat);
					}
					else {
						MPI_Recv(&x[it][(i-1)*tam],tam,MPI_DOUBLE,i,1,MPI_COMM_WORLD,&stat);
					}
				} //recebe o x calculado por cada filho

				//acha o maior valor de Xk+1 - Xk (1 por 0)
				max1 = fabs(x[it][1] - x[it-1][1]);
				max2 = fabs(x[it][1]);
				
				for(i = 1; i < J_ORDER; i++) {
					if(max1 < fabs (x[it][i] - x[it-1][i]) ) {
				    	max1 = fabs (x[it][i] - x[it-1][i]);
					}

					if(max2 < x[it][i] ) {
				    	max2 = fabs(x[it][i]);
					}

				}

				e = max1/max2;
				if(e <= J_ERROR) {  //se for menor o erro, continuar recebe 0 e sai do while 
					continuar = 0;
				}
		    	
				for(i = 1; i < numTasks; i++) {
					MPI_Send(&continuar, 1, MPI_INT, i, 1, MPI_COMM_WORLD); //envia para os filhos o sinal de continuar ou não.
					MPI_Send(&x, J_ITE_MAX*J_ORDER, MPI_DOUBLE, i, 1, MPI_COMM_WORLD);
					if(i == numTasks-1) {
						MPI_Send(&tamUltimo, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
					}
					else {
						MPI_Send(&tam, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
					}
				}
			}
			
			for(i = 1; i < numTasks;i++){
					MPI_Recv(&it, 1, MPI_INT, i, 1, MPI_COMM_WORLD, &stat);
			}
	
			comprovacao(J_ORDER, x[it], J_ROW_TEST, matriz, b, it);
			
			fim  = clock();
			double  tempo = (fim - inicio)/(CLOCKS_PER_SEC/1000);
			printf("tempo = %g\n", tempo);
		   
		}
		else {
		    int continuarFilho, itFilho, tamFilho;
		    double xFilho[J_ITE_MAX][J_ORDER];
		    
		    //recebendo as informacoes do pai
		    MPI_Recv(&xFilho, J_ITE_MAX*J_ORDER, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD, &stat);
		    MPI_Recv(&continuarFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &stat);
		    MPI_Recv(&itFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &stat);
		    MPI_Recv(&tamFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &stat);

		    while(continuarFilho == 1) {
				itFilho++;
				
				for(i = (myRank-1)*tamFilho; i < myRank*tamFilho; i++) {  
					for(j = 0; j < J_ORDER; j++) {
				    	if(i != j) {
				        	xFilho[itFilho][i] += (aEstrela[i][j] * xFilho[itFilho-1][j]);
						}	 
					}
				    xFilho[itFilho][i] = xFilho[itFilho][i] + bEstrela[i];
				}
				
				MPI_Send(&itFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);
				//envia para o pai tudo o que foi calculado
				MPI_Send(&xFilho[itFilho][(myRank-1)*tamFilho], tamFilho, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
				//espera o sinal de continuar
				MPI_Recv(&continuarFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &stat);
				MPI_Recv(&xFilho, J_ITE_MAX*J_ORDER, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD, &stat);
				MPI_Recv(&tamFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &stat);
			}
			MPI_Send(&itFilho, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);

		}
    }
    
    #pragma omp parallel for
    for(i = 0; i < J_ORDER; i++) {
    	free(matriz[i]);
    }
    free(matriz);
    
    #pragma omp parallel for
    for(i = 0; i < J_ORDER; i++) {
       	free(aEstrela[i]);
    }
    free(aEstrela);

	fclose (pFile);
    
    MPI_Finalize();
		
	return 0;

}
